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ABSTRACT 

We obtained self-similar solutions of relativistically expanding magnetic loops by as- 
suming axisymmetry and a purely radial flow. The stellar rotation and the mag- 
netic fields in the ambient plasma are neglected. We include the Newtonian gravity 
of the central star. These solutions are extended from those in our previous work 
( Takahashi. Asano. fc Matsumoto 2009f ) by taking into account discontinuities such 



as the contact discontinuity and the shock. The global plasma flow consists of three 
regions, the outflowing region, the post shocked region, and the ambient plasma. They 
are divided by two discontinuities. The solutions are characterized by the radial ve- 
locity, which plays a role of the self-similar parameter in our solutions. The shock 
Lorentz factor gradually increases with radius. It can be approximately represented 
by the power of radius with the power law index of 0.25. 

We also carried out magnetohydrodynamic (MHD) simulations of the evolution 
of magnetic loops to study the stability and the generality of our analytical solutions. 
We used the analytical solutions as the initial condition and the inner boundary con- 
ditions. We confirmed that our solutions are stable over the simulation time and that 
numerical results nicely recover the analytical solutions. We then carried out numer- 
ical simulations to study the generality of our solutions by changing the power law 
index 5 of the ambient plasma density po r^^ . We alter the power law index 5 from 
5 ~ 3.5 in the analytical solutions. The analytical solutions are used as the initial con- 
ditions inside the shock in all simulations. We observed that the shock Lorentz factor 
increases with time when the power law index is larger than 3, while it decreases with 
time when the power law index is smaller than 3. The shock Lorentz factor Fs can 
be expressed as Ts oc t^^~^'^/'^ where 5 is the power law in dex of th e amb ient plasma. 
These results are consistent with the analytical studies bv IShapirol ( 19791 ). 
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INTRODUCTION 



Soft Gamma-ray Repeaters (SGRs) and Anomalous X-ray Pulsars (AXPs) a re believed to be youtiR neutron stars with str ong 



Woods fc Thompsonll2006l : iMereghettil boosi , for 



magnetic fields (~ 10 G). They are categorized as magnetars (see, e.g., 
review). Their rotation period P and its time derivative P are P ~ 10 s and P ~ 10~^° s s~^, respectively. A magnetic 
field strength inferred by assuming the dipole emission from P and P is about 10^^ Gauss. Persistent X-ray emissions with 
the luminosity of Lx ~ 10"^* — 10^^ erg s~^ are observed in SGRs and AXPs. SGRs are identified by the hard X-ray bursts. 
Extraordinary energetic outbursts called giant fiares are observed in three SGRs. The burst energy in the SGR 1806-20 giant 
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flares on December 27, 2004, reached ~ 10*^ erg ( Terasawa et al.ll2005f ). Since SGRs and AXPs are not accretion powered 
sources, it is befleved that their activity is driven by dissipation of their magnetic fields. 

Strong magnetic fields are created by the dynamo mechanism during the core collapse of a supernova progenitor. At this 
stage, the star becomes unstable again s t the convective mode since the entropy gradient becomes negative, (dS/dr < 0) due 
to the neutrino cooling (Burrow^ 1987 : Keil et al. 1996h . Thus the infant neutron star can store a large amount of magnetic 
energies. The magnetic fields are also amplified after the birth of the magnetar. The interior of the star rotates differentially 
at its birth when the equation of state is stiff. The internal magnetic fields can be amplified up to lO^^G due to the dynamo 
mechanisms IjDuncan fc Thompsonlll992r ). As the magnetic helicity is accumulated inside the star, the Lorentz force exerted 
by the twisted magnetic fields balances with the rigidity of the crusts. When the critical twists are accumulated, the crustal 
rigidity can no longer sustain the Lorentz force by the strong magnetic fields. The magnetic helicity is then injected into 
the magnetosphere. The resulting crustal motion induces the electric fields and it results in creating the potential difference 
between the foot-points of the magnetic loops. The particles injected from the interior of the star are accelerated along the 
magnetic field lines due to the p otential difference. The accelerated particles initiate the avalanches of the pair creation 
i Beloborodov &l ThompsonI 120071 ). These particles carry the electric current, which twists the global magnetic fields. When 
the dynamical equilibrium is lost by the accumulated magnetic twists, the magnetic loops expand by the magnetic pressure 
gradient force. Inside the magnetic loops, a current sheet similar to that of solar flares is formed. The magnetic recon nections 
taking place in the current sheet are responsible for the magnetic energy release and resulting flares ( Lvutikov 20061 ) . Recent 
observat ions which indicate the topological change of the global magnetic fields before and after the giant fiares support these 
models (jWoods et ahlboOlh . 

Motivated by the magnetar flare model, Spitkovskv 1 20051 ) performed 2-dimensional relativistic force- free simulations of 
magnetar flares by injecting the magnetic twists into the magnetosphere. They showed that the initially dipole magnetic fields 
are tw i sted b y the foot-point motion and the loop magnetic fields then expand due to the magnetic pressure gradient force. 
Asand l|2007l ) carried out 2-dimensional relativistic force-free simulations of exp anding magne tic loops and showed that the 
Lorentz factor defined by the drift velocity Vd = c{E x B)/B^ exceeds 10 (see. IUchidalll997l . for the definition of the drift 
velocity). These simulations indicate that the magnetic loops expand self-similarly. 

S uch self-similar solutions have been found in analytical studies. In the framew ork of the forc e -free d ynamics, iLvutikov fc BlandfordI 
( 2003h obtained self-similar solutions of the spherically e xpanding magnetic shell. |P rcndcrgast| (|2005n found self-similar solu- 
tions of the relativistic force-free field in two dimensions. Gourgouliatos fc Lvnden-Bell (,200^ derived relativistic s elf-similar 
force- free solutions and analyzed them in detail. In the framework of the relativistic magneto hydrodynamics fMHD) , iLvutikov 
I 2OO2I ) bund self-similar solutions of the spherically expanding magnetic shells. Recently, I Takahashi et al.l (120091) obtained 
self-similar solutions of magneti c loops (not shell) by extending the non-relativistic solutions obtained by Low ( 1982lJ ). 
Gourgouliatos fc VlahakisI (|2O10l ) obtained solutions by ignoring the gravity from the central star. These authors studied 
outflows of the the magnetized plasma hfted up from the central star. However, they did not consider the interaction between 
the outflow and the interstellar matter. iLowl l|l984af ) obtained non-relativistic self-similar solutions of the expanding magnetic 
loops interacting with the interstellar matter. In their models, the outflows and the ambient plasma are divided by a contact 
discontinuity. The forward propagati ng wave forms anot her discontinuity (shock). This solution is useful to understand the 
coronal mass ejections in solar flares. Stone et al. 1 1992h employed this solution as a test problem to check the validity and 
accuracy of axisymmetric MHD codes. 

Takahashi et al.l (|2009l ) by including the contact discontinuity 



In this paper, we extend the analytic solutions give n by 



and the shock by extending the non-relativistic model by [Low 



1 1984al ) to the relativistic regime. 



This paper is organiz e d as follows. In § [2l we summarize the basic equations of self-similar relativistic MHD equations 
given bv I Taka hashi et"aL (|2009|). In § |3] we show the solutions of these equations including the two discontinuities. These 
solutions represent the relativistic coronal mass ejection from the central star. Such solutions are expected to explain the 
giant flares in magnetars. The physical properties of the solutions are shown in § [l] We also carried out the 2-dimensional 
relativistic MHD simulations to study the stability of the solutions. The analytical solutions shown in § [3] are used as the 
initial and boundary conditions for simulations. These results are shown in § [S] We summarize our results in § [S] 



2 BASIC EQUATIONS OF SELF-SIMILAR RELATIVISTIC MHD 



In the following, we take the light speed as unity. The complete set of the relativistic ideal MHD equations is 
d 



dt 



{ip) + V ■ (7pv) = 0, 



PI 



4- (w ■ V) 



dt 



+ (w • V) 
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V - B = 0, 

at 



47rj, 



(4) 
(5) 

(6) 

(7) 
(8) 



where E, B , j jV,'}, pe, p,p,T are the electric field, the magnetic field, the current density, the velocity, the Lorentz factor, the 
charge density, the mass density, the pressure and the specific heat ratio, respectively. The vector er is a unit vector in the 
radial direction. Newtonian gravity of the central star is included as an external force. Here G is t he gravitational constan t, and 
r is the distance from the centre of the star. We replaced the gravitational force ~GM p-y /r^ in Takahas hi et al.| ( 200E ) with 
—GM phj^/r^, which can treat the gravitational force of the relativistic plasma more properlv ijCourgouliatos fc VlahakisI 
2O10l ). 



The relativistic specific enthalpy including the rest mass energy h is written as 



h 



1 + 



P 



1 + /ijv. 



(9) 



p r-ip 

where e is the energy density of the matter including the photon energy coupled with the plasma, and Hn ~ rp/[(r — l)p] is 
the non-relativistic specific thermal enthalpy. 

In the following, we take F = 4/3 which corresponds to the relativistic radiation pressure dominant plasma. Thus we can 
study the evolution of a fireball confined by magnetic fields. 

We ignore the stellar rotation and assume a purely radial flow. We also assume ajdsymmetry. The magnetic fields in 
spherical coordinates (r, 9, (j>) are then expressed in terms of two independent scalar functions A and B as 



1 dA dA 



B 



(10) 



r sin 9 \r 86 ' dr 

In the following, we assume that the evolution of the magnetic loops can be described by a Lagrangian coordinate 77: 

(11) 



r 

wy 



where Z(t) is a scale function of time. We further assume that the flux function A evolves with time t and radial distance r 
through the Lagrangian coordinate r), as 



A{t,r,9)^ A(ri,6). 

The MHD equations are then written as 

dt2 

Vr = Z{t)ri, 
pit,r,9) = 

p(t,r,0)7 
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is an operator: 
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(see, iTakahashi et alJboogl . for derivation). We can readily solve equation (|13|l as 



m = ^^t, (22) 

where ^ is a constant of integral. The radial velocity is obtained from equation (|14|l as 

Vr = VIti- (23) 
By substituting equation (|22l) into equation (|19p . we obtain the function Q which has a form: 



«('''^) = rr^' (24) 

where (?(A) is an arbitrary function of A. 

We note that the last term in equation (|18|) includes Z{t) = y/^t, which is a function of time. Thus this term does 
not admit the existence of the self-similar solutions because the thermal enthalpy in the gravity term explicitly depends on 
time. This violation of the self -similarity comes from the differ ence of the scaling law between the density and the pressure 
(equations 1151 and 1161 also, see Gourgouliatos fc Vlahaki^ 2O10l ). However, when the contribution of the thermal enthalpy on 



the gravity is sufficiently small, we can obtain self-similar solutions of the relativistic expansion. The ratio of the gravity for 
the thermal enthalpy to the plasma inertia is written as 

gravity for the thermal enthalpy _ \GM phNj'^ /r^\ _ rg 
plasma inertia \p'yD{'yhv) / Dt\ ru^ ' 

where rg — GM is the gravitational radius. Thus when r ^ rg, the gravitational force for the thermal enthalpy /ijv is negligible 
when « ~ 1. In this paper, we consider the evolution of the relativistically expanding plasma at large distance where the 
gravity for the thermal enthalpy —GMphN^^/r'^ can be neglected. In such region, the plasma inertia is sustained by the 
pressure gradient and Lorentz forces. When r ^ _Rs, we assume that the rest mass energy density much exceeds the thermal 
energy density (i.e., ht^ <^ 1). In this regime, the gravity is expressed as —GMp')^/r^. The Lagrangian coordinate r] (see 
equation lll[) then behaves as the self-similar parameter. 

Without loss of generality, we can change independent variables of P from (r), 9) to (r), A) as 

P(77,e) = P(r,,i). (26) 

By neglecting the thermal enthalpy in the gravitational force —GMph'y^/r'^, equations (|20p and (|18l) reduce to the self-similar 
equations as 



dP _ 1 
dA ^irri'^ sin^ 8 

D 



p , d (, 2dA\ , 9iA) dg{A) 



dri \ drj J 1 - ^r]^ dA 



rf\P^^iff ( 4:^vP dP 



GM dr) 



(27) 
(28) 

Equations (|24l) . (|27|) and (|28[) are the set of the self-similar MHD equations. The explicit solutions can be constructed as 
follows. First we prescribe an arbitrary function A{r],6) and the function g{A). The pressure function P is determined from 
equation (I27p and the density function D is obtained from equation (|28p . 

Before presenting the solutions of the self-similar MHD equations, we have to note that the self-similar relativistic ideal 
MHD equations describe the free expansion of the magnetized plasma. By taking time derivative of i), we obtain 

where we used equations (llip . (|22|l and (|23[l . The plasma is neither accelerated nor decelerated, but it expands with the 
inertial speed keeping the for ce balance (this can be confirmed by inserting equation (I23|) into equation ((2]), or, see equation 
(36) in ITakahashi et akllioogl ) . 



3 RELATIVISTIC CORONAL MASS EJECTION 



In the previous section, we showed the set of the self-similar relativistic MHD equations, which describes the plasma expanding 
with the inertial speed. In this section, we obtain solutions of these equations by imposing a ppropriate boundary conditions. 

We adopt the simple model of magnetic explosions according to Low models (|Lowlll984al . see, Fig.[l|. Equations (I27|l and 
(|28|l describe magnetized plasma outflowing from the central star (outflow region). The outflow sweeps up an ambient plasma 
while it expands. A contact surface would be formed which separates the outflowing plasma with the swept-up ambient plasma. 
The contact surface is situated at r = Rc(t). Ahead of the contact surface, a forward wave propagating into an undisturbed 
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Figure 1. Schematic picture of our models of the relativistic coronal mass ejection. Ahead of the global magnetic loops, two discontinuities 
exist at r = Rc{t) (contact discontinuity) and r = Rs{t) (shock). 



ambient plasma forms a shock at r = Ra{t). The swept-up ambient plasma is accumulated in the region in Rc{t) ^ r ^ Rs{t) 
(postshock region). 

First, we consider the outflow region where the magnetized p lasma i s lifted up from the c entral star. Various models of 
the magnetic field confi gurations ha ve been proposed (Low..l982bl . Il984al ; iTakahashi et al.lboogi ') . In this paper, we utilize the 
following fiux function l|Lowlll984al ) 



Ao 



where 



Aof{Xr]) sin 



fix) ^Ho + 



(30) 



(31) 



Here Ao, Hq, and A are constants. The fiux function given by equation pop has local maxima (see Fig. 4 in lLowlll984af ). The 
maximum corresponds to the centre of the fiux ropes. This function can thus describe the cor onal mass eject ion from the 
central star. Such relativistic coronal mass ejections can be applied to giant fiares in magnetars ( Lvutikov 20061 ). We express 
the toroidal magnetic field gr as a power series of A: 



5(i) = ^a„l", 



(32) 



where are constants. Note that the dependence of the toroidal magnetic field on the polar angle is represented by 
sin^"~^ S from equations (|10|) . (|17|) . (|30|l and (|32}. This indicates that the power law index n and a„ can be inter- 
preted as the Fourier modes of the toroidal magnetic fields in polar direction and their amplitudes, respectively. The modes 
and the amplitudes are determined by the boundary condition at which the magnetic shear is injected into the magnetosphere. 
The shear injection begins before the self-similar expansion starts. In this paper, we do not consider the details of the shear 
injection from the central star and leave them as free parameters since we consider the self-similar stage. 
The explicit forms of the magnetic fields are then written as 



sin(A77) 



— cos(A?7) 



cos 8, 



(33) 



Bg{t,r, 



4V- 

1 



cos(A77) — ^ — ^^^^ sin(Ar7) 



Ar? 



smfc'. 



(34) 
(35) 



sin 9 1 — ^rj^ 

n 

The pressure function of the outfiowing plasma Po is obtained by inserting equations (|30p and (|32p into equation (|27p as 
Po{r,,e) ^ PA{il,e) + PQ{ri,e) + P(ri), (36) 



where 
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Pa{v,0) = 



2Ho + J- (Ar,) [(1 + Cr?') sm(Ar,) - {\rj){l - C^?') cos(A,,)] 



(37) 



E 



Anf]'^ ( 1 — ,f 7/2 ) I m + n — 1 



(Ar/jsm 



2{m+n-l) 



'+ ^ mamariAof{Xri)lnA{-q,9) 

m + n— 1 J 



(38) 



Here Pi arises from the integration. By substituting equations (|36 
Do is expressed as 

Doiv,e) ^ DA{v,e) + DQ{r,,e) + D,{r,), 
AoA{Tj,9) 



into equation 



DAiv,0) 



--^iv), 



the corresponding density function 

(39) 
(40) 



DQ{v,e) = — 



i/o + Vf(^'/)sm(Ar?) 



E ■ 



m + n — 1 



47rGMr7(l -C»72)2 



mamQ:„Ao In A(r;, 61), 



for m + n 7^ 1, 



(41) 



m+n — 1 



for m + n = 1, 



where the function 'S(rj) is defined as 



Sin) = 8//o + W -Ar, { [(3 - (A,,)^)(l + (C^?')') + 2Cr?'(l + (A,,)^)] sin(A7;) - A77(l - ^77^)(3 + 5^v^) cos{\rj)} 



(42) 



The function Di{rf) describes the isotropic distribution of the plasma which is related with Pi{ri) through equation (|28|) . 

To determine functions Pi and Di, we need one more relation between them. We assume that there is no energy/mass 
injection from the central star at this stage. Since the contact surface separates the outflow region from the postshock 
region, the outflow plasma expands adiabatically. We then obtain another relation between Pi{ii) and Di(rf) from the entropy 
conservation equation as 

^ = const (43) 

where pi — Z{t)'^Pi, pi = Z{t)'^Di and is a constant. Substituting equations psp . (jlGp . and 1)23^ into equation (|43|l . we 
obtain 



P 



2^ 3 



(44) 



Substituting equation (|44l) into (|28[) . we obtain the solutions. 



(45) 



A(r?) = 



(46) 



where ^ is a constant of integral. These solutions describe the isotropic outflowing plasma in the outflow region. 

Next, we consider the post shock region {Rc = r ^ Rs)- The shocked plasma also moves in radial direction. We assume 
that the shocked plasma evolves self-similarly and obeys the same set of self-similar equations in the outflow region. Then the 
contact surface moves in radial direction with a constant speed and the radius of the contact surface Rc is expressed as 



(47) 



where 77c is a constant. We assumed that the magnetic fields in the ambient plasma can be neglected. From these assumptions, 
the shocked plasma obeys equation (|28[) . We need another relation between the shocked gas pressure Ps and the shocked gas 
density Ds- Note that we cannot use the adiabatic relation because the ambient plasma flows into the postshock region from 
the shock surface at r = Rs. The forward shock compresses and heats up the plasma, resulting in an increase in the entropy. 
Thus the entropy in the shocked plasma should be determined by the shock condition at r = Rs{t). Rather than evaluating 
the entropy variation by the shock, we consider the jump conditions of the plasma density, the pressure, and the velocity. The 
entropy vari ation is determined after imposing the Rankin-Hugoniot relations between the undisturbed and shocked plasma 
I Low 



py van ; 
ll984bll 
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By assuming the strong shocks, the relativistic Rankin-Hugoniot relations are written as 



2p2 



P'y\r = R 
2 I 



2r^/9o| 



7 



r=Rs 2 



(48) 
(49) 
(50) 



with an accuracy of 0(1/72) (jBlandford fc McKee|[l976l : iKennel fc Coronitilll984r ). Fs and po are the shock Lorentz factor 
and the plasma mass density of the undisturbed ambient plasma, respectively. Here we ignore the thermal pressure in the 
undisturbed plasma by assuming the strong shock. 

From equations (|23p and (|50p . we obtain the time evolution of the shock radius as 

SRoXs 



Rs{t) = ■ 

where Xs 
and Ds'. 



= Rs/t and Ro is a constant of integral. Combining equations 



Here we use equations p5|) . (|16|) . (|51|) . Substituting equation (|52p into equation (|28p . Ps and Dj, are expressed as 
Po 



(51) 

and (|49|) . we obtain the relation between the Ps 

(52) 



Ps = 



V2 



where Po is a constant of integral. The function K is expressed as 



K 



; exp 



GM 

' Ro 



dx- 



+ VTT^) (1 - : 



Po{r 



(^/2x + VTT^) 

The ambient plasma density is obtained by substituting equations (|51|) and (|54p into (|49p as 
3Po Vo 



-K 



(53) 



(54) 



(55) 



(56) 



(57) 



4r4 (1 - ^r,l) 

Here rjo should be determined from the following equation, 
3Ro\^rjo 

r — 

[V^VO + VTT^) (^/2Cr?o - yr+C^) 

Note that the radial profile of the ambient plasma density is not arbitrary but determined by equation (|56p . Some authors 
derive the self-similar solutions by p rescribing the density profile of the ambient plasma as po oc r"*, where 5 is a constant 
( Blandford fc McKedl 19761 : [Sajill2006 h . ] n our approach, we first prescribe the self-similar variables given in equation (jlip . Then 
the outflow velocity is obtained by equation (|14p . The ambient plasma density is determined by applying the Rankin-Hugoniot 
relations at the shock. Thus the ambient plasma density cannot have an arbitrary form, but it is uniquely determined. 

Finally, we apply the boundary conditions at r = Pc. Since we assumed unmagnetized ambient plasma, the magnetic 
fields should vanish at r = Pc. This condition can determine the parameters Hq and A in equation (|30p . The conditions that 
Br{r — Rc) = and ^^(r — Rc) = are expressed as 



P^o 



sin(A7;c 



cos(A77c 



from equations (|33l) and psp . Another condition is that Be{r — Rc) — 0. This condition is written as 



tan(A77) = 



(Ar,)2 



(58) 



(59) 



Equation (|59|l determines A and then the parame ter Po is obtained from equation (|58|) . Note that the equation (|59p has an 
infinite number of roots (see Fig. 4 in lLowlll984al ). The first root for rj > arises at Xrj — Xrji ~ 2.7 and the second root does 
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Figure 2. Contour plots of the magnetic flux A (left) and the toroidal magnetic fields. The contact surface is situated at rj/rjc = 1. The 
parameters are taken as r)c = A, ^ = 1 and n = 1. 



at Xi] = \ri2 ~ 6.1. The first root correspon ds to the po sition of the centre of the flux ropes. We take the second root as the 
contact surface, i.e., rjc = ')2 {zi in Fig. 4 of |Low|[l984a|) throughout this paper. The parameter Ho is then determined from 
equation (|58|l as Ho — 0.81. 

Another constraint on the parameter comes from pressure balance across the contact discontinuity. The pressure Po 
consists of three component of the pressure. Pa, Pq and Pi. Pa and Pq are exactly zero ai r — Rc since A\r,^,^^ = dA/dr]\n^r]^ = 
. Pi should be smoothly connected with Ps at the contact surface. From this condition, the parameter Po is expressed as 



(60) 



Here we used equations (|45|l and (|53[) . 

The remaining parameters are r]c, Ro, otn, n, v, and /i, where ^ denotes the scaling of time and radius and 77c describes 
the velocity of the contact surface. Equation (|51|) determines _Ro by prescribing the shock radius when self-similar expansion 
starts. The twist injection at the central star determines the amplitude an and the Fourier mode number n of the toroidal 
magnetic fields. A constant v which appears in equation 1)43^ denotes the entropy of the isotropic plasma in the outflow region 
r ^ Re- The density at the contact surface when self-similar expansion starts determines the constant fi which appears in 
equation (|46}. 



4 PHYSICAL PROPERTIES OF THE SELF-SIMILAR EXPLOSIONS 

First, we concentrate on the structure of the magnetic loops in r ^ Rc. Fig. [2] shows contours of the magnetic flux A (left) 
and the toroidal magnetic fields (right) in the i] — 6 plane. The parameters are rjc = X and ^ = 1. The poloidal mode number 
of the toroidal magnetic fields is n = 1. We can see the fiux rope structures emerging inside the expanding magnetic loops. 
The centre of the flux ropes is situated at rj — rji (771 ~ 0.4477c). Since the magnetic fields have both poloidal and toroidal 
components, they describe the twisted fiux ropes. The fiux ropes rise in the +r direction with time. In the limit of t 2> r, the 
magnetic fields are represented as 

lim B — ^"^P^'^ cos 9 , (61) 
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Figure 3. Contour plots of the pressure (left) and the density (right). The contribution from the isotropic plasma is subtracted. The 
contact surface is situated at ri/ric = 1. The parameters are taken as Jjc = A, ^ = 1, n = 1, and fi = 0. 



from equations (|33|) - H35[) . The toroidal magnetic fluxes are diluted by the expansion according to the flux conservation 
equation. The configuration of the magnetic fields approaches that of the split monopole. 

Fig. [3] shows contours of the gas pressure (left) and the density (right) in the rj — 6 plane. Contributions from the 
isotropic plasma, Pi and Di, are subtracted, so that the pressure and the density can be negative in these panels. The set of 
the parameters is the same as those in Fig. [2] Accompanying the flux ropes, low density voids are generated. The toroidal 
magnetic fields can create such voids. As mentioned at the end of § [2] since the force balance is attained in our self-similar 
solutions, the Lorentz force by the poloidal magnetic fields should balance with the gas pressure gradient force and the toroidal 
magnetic pressure gradient force. This indicates that the gas pressure decreases as the toroidal magnetic pressure increases. 
Such voids exert the buoyancy force on the plasma in radial direction. Fig. [3] shows that the pressure gradient force ahead of 
the voids balances the buoyancy force. 

Next we consider the isotropic part of the outflowing plasma expressed in equations (|45|) and (|46p . When the constant 
fx is exactly zero, the plasma distribution reduces to that in hydrostatic states. The scaling comes from the assumption of 
adiabatic expansions with the polytropic index of F = 4/3. When ^ is not zero, the plasma distribution differs from that in 
the hydrostatic st ates. When t he flow speed is non-relativistic, i.e., y/S^rj <^ 1, the solutions reduce to those in non-relativistic 
MHD obtained by Low ( 1984a ). Since the enthalpy contributes to the plasma inertia, the correction term (1 — ^rj^)~^^^ arises 
in the relativistic MHD. The plasma tends to be hydrostatic since ^/^rj <^ 1 when t^r. 

In front of the outflow region, the ambient plasma is compressed by the shocks at r = 7?s and is accumulated in the post 
shock region {Rc ^ r ^ Rs). The contact surface at Rc = \^'nct divides the outflowing plasma from the shocked plasma. 
Since we assumed that the postshock gas evolves self-similarly according to the same basic equations for the outflow region, 
the contact surface has the constant velocity Vc — \/^r]c. The expansion speed of the shock radius dRs/dt is, however, not 
constant but it increases with time. The time evolution of the shock radius is expressed by equation (|5ip . Fig. |4] shows the 
radial profile of the outflow Lorentz factor. Dashed, dotted, and dot-dashed curves denote the Lorentz factor ai t — 10, 25, 40, 
respectively. Thick curve denotes the time evolution of the shock Lorentz factor Fs, while thin solid curve does the time profile 
of the outflow Lorentz factor at the shocks (jlr^R^)- The value of the parameter Ro is Ro = 1.78 x 10"*, which corresponds 
to 'y{t = 5) = 8. The shock Lorentz factor is larger than the outflow Lorentz factor at r = Rs by factor \/2, as expected 
from relativistic strong gas dynamical shocks (see equation I50p . The undisturbed plasma is abruptly heated up by the shocks. 
Plasma velocity suddenly becomes zero ahead of the shocks (r > Rs) where the undisturbed plasma exists. 

As the shock propagates in the undisturbed plasma, the shock surface is accelerated. The shock Lorentz factor and the 
shock radius can be expressed as 
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Figure 4. Radial profile of the outflow Lorentz factor. Dashed, 
dotted, and dot-dashed curves show the outflow Lorentz factor 
aX t = f 0, 25, 40, respectively. Thick solid curve shows the shock 
Lorentz factor, while thin solid curve does the maximum Lorentz 
factor of the outflow. The value of the parameter i?o is taken to be 
Rq = 1.78 X 10^^, which corresponds to 7(t = 5) = 8. 
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Figure 5. Radial proflles of the plasma pressure (thick solid curve), 
the plasma density (thin solid curve) and the Lorentz factor (dashed 
curve). Radius is normalized by the shock radius Rs- The contact 
discontinuity is situated at Rc/Rs ~ 0.834. We take Ts = 8\/2 and 
/i = 0. 
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respectively. Here we approximate w ~ 1 to obtain the first equality in equation (|62p . While the shock surface moves with 
almost constant speed ~ 1, the shock Lorentz factor increases with time with the power law index of 0.25. This result comes 
from the density profile of the undistu rbed plasma. As shown later, the undisturbed plasma density decreases with radius as 



r~^ ''^ (see equation [64]) . Shatairo ( 19791 ) showed that when p cc r"", the fiow is acc elerated when S > 3. Although we take into 
account the gravity from the central star, which is not included in Shapirol ( 19791 ). the gravitational force is smaller than the 
pressure gradient force in this region. Thus their relation can be adopted in our analysis. The outflow is gradually accelerated 
when 5 ~ 3.5. 

Next we consider the radial profile of the plasma pressure and the density. Fig. [5] shows the radial profile of the pressure 
p (thick solid curve), the density p (thin solid curve) and the outfiow Lorentz factor (dashed curve) on the equatorial plane. 
The horizontal axis denotes the radius normalized by the shock radius R^. The shock Lorentz factor is taken as — 8\/2. 
We take ^ = for simplicity. The other parameters are ^ — I, Ro — 1.3 x 10"*, Qn = 0, and Ao = 0. The contact surface is 
situated at Rc — 0.834_R3. Behind the contact surface, the gas pressure and the density decrease with radius as oc and 
oc r"'^, respectively. The density jump appears at r = Rc, while the pressure is continuous (contact surface). The swept-up 
ambient plasma is accumulated in region Rc ^ r ^ Rs- The plasma density in the postshock region (0.834 < r/Rs < 1) is 
larger than that in the outflow region. The strong pressure gradient force in this postshock region pushes the plasma in +r 
direction. The pressure gradient force balances with the inertia of the accumulated plasma. Thus, the plasma flows toward 
+r direction with the inertial velocity. 

Ahead the contact discontinuity, a strong shock appears a,t r — Rs- We assumed the strong shock and neglected the 
ambient plasma pressure. The plasma is abruptly heated up by the shock. The ambient plasma density jump also appears at 
the shocks. The density of the shocked gas is larger than that of the ambient plasma by factor 32 for Fs = 8\/2 as expected 
from the relativistic Rankin-Hugoniot relations (see equations 1491 and I50() . The density of the undisturbed plasma is described 
by equation (|56p . From this equation, the density can be approximately represented by the power law of r as 



p{r) 
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(l + ^/2) 
6 



V2 



-A'(l)r" 



(64) 



Here we approximate Vr{r — Rs) — 1 for relativistic flows. The density decreases with radius with the power law index of 



3.5. The ambient plasma density decre ases slightly faster than that inside the shocks (oc 



The shock Lorentz factor 



thus increases with radius (|Shapirolll979l ). 

Next we consider the total energy contained within the spherical surface Rs . Let £ be the total energy in r < Rs 
shown in Takahashi et al. I l|2009l ). the virial theorem can be applied for the relativistic inertial flow: 



As 
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£ = K. ~ (Sr - 4)Uti, + J ^^j;^dV + J pr-dA--^ J {2 [{r ■ E){E ■ d A) + {r ■ B){B ■ d A)] - {E^ + B^){r ■ d A)} , (65) 
where 

IC= dVp7^ (66) 



f/th = J dV^, (67) 

and S shows the Poynting flux. Here A denotes expanding spherical surface at r = Rait). The third, fourth, and fifth terms 
of the right hand side of equation (|65p represent the Poynting fiux, the work done by the gas pressure, the work done by the 
Maxwell stress, respectively. 

Let us evaluate the non-kinetic part of the energy, £' — £ ~ IC. The second term in the right hand side of equation (1651) 
is zero because the polytropic gas index is F = 4/3. The fifth term of (I65() is zero since the electromagnetic fields vanish at 
r = Ra. The third term also becomes zero after the straightforward calculations. Thus the non-kinetic part of the energy can 
be evaluated as 



£'{t) = ^poT'iRl 



47r „Q 



(68) 



Here we used equations (|48p and (|50[l . The energy £' does not depend on the amplitudes of the magnetic fields. As shown in 
§[31 the plasma density (pressure) consists of three parts. Da, Dq and Di {Pa, Pq and Pi). While the components Da and 
Dq depend on the magnetic field strength, Di is independent of them. Da and Dq do not contribute to the non-kinetic part 
of the energy from equation (|65|) since they are exactly z ero at r = _Rc- This indicates that the plasma interacting with the 



magnetic fields is in marginally stable state. iLowl (|l982al ) showed that the inertial flow with the polytropic index F = 4/3 
represents the marginally stable state in non-relativistic MHD. This situation is also valid for the relativistic MHD. Total 
energy contained inside the shocks is thus independent of the strength of the magnetic flelds. 

By substituting equations (|62ll and (|64p into equation (|68p . the time evolution of the energy £' is written as 

S'^(l±f}Zm^, (69) 
6 4 -Ko 

where we used equations (|62p . (|63p . and (|64[) . Note that the non-kinetic part of the energy £' is positive. This means that 
plasma speed exceeds the escape velocity determined by the gravitational potential. 
The kinetic energy IC given in equation (|66l) is expressed as 



IC^2tv I desine I dri I^^S^^ (70) 



where we used equations (|16|1 and (|23|) . Note that equation (|70|l depend on time through 77. When we integrate inside the 
sphere with radius r = Rs, the integration is carried out in [0, r^s] in the self-similar space. Since the flow is relativistic, i.e., 
Vr ~ 1, r^s — l/\/C is constant with time. Thus both non-kinetic and kinetic energies are constant with time. Strictly speaking, 
the total energy should increase with time since the shock surface sweeps up the ambient plasma. The rest mass energy of 
the swept-up plasma contributes to the increase in the total energy. This energy is, however, negligible because we assume 
the strong shocks. From equation (|48() . the rest mass energy density of the undisturbed ambient plasma is smaller than the 
kinetic energy density of the shocked gas by factor F* and negligibly small (Note that the relations [48] - [50] are correct with 
the accuracy of 0(1/F^)). For the same reason, £' is also independent of time with the accuracy of ©(l/F^). Using these 
facts, the shock Lorentz factor is expressed as 

l.£'^p-^[r = Ra)R:^, (71) 



from equation (|68|) . This result is equivalent with equation (16) in Shapirol 1 19791 ). Although our solutions include the magnetic 



flelds and the gravity from the central star, they do not contribute to the total energy because the system is in a marginally 
stable state. Thus only the hydrodynamical (isotropic) part contributes to the the total energy. By inserting equations (|63[l 
and (|64l) into this equation, we obtain Fs oc r^^'^ again (see equation I62[). The shock is thus accelerated when it propagates in 
the ambient plasma. 



5 NUMERICAL SIMULATIONS 

In this section, we show results of relativistic MHD simulations to study the stability of our solutions. For this purpose, we 
use the analytical solutions as the initial conditions of the numerical calculations. The relativistic MHD equations are solved 
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Figure 6. Colour contour shows the gas pressure profile and curves show the magnetic field lines for analytical solutions (left) and 
numerical results (right). 



in two dimensions using polar coordinates (r, 6). We assume axisymmetry {d/{d4>) ~ 0). The number of grid points is (A'^,., 
Ne) — {3600, 360) on the domain of R^n = 1 ^ r < 50 in normalize d unit and ^ 6 ^ tt. The grid sizes are Ar = 1.39x10"^ and 
Ad = 8. 72 X 10~^. We use the H LL method ( Harten et al. 19831 ) to calculate numerical fluxes. We utilize the modified CTU 
method I Mignone fc Bodol 20061 ) to achieve the second order accuracy in space. In our analytical solutions, the strong shock 
is expected. Such a strong shock can indu ce the numerical oscillations. To avoid the problem, we utilize the harmonic mean 
for smoother prescription (jvan Leeii 119771 ). We use the Constraint-Transport method to satisfy the no monopole condition. 
We impose the outflow boundary conditions in outer radial boundary at r = 50 and the symmetric condition at the axes at 
6 — 0,n. The inner boundary conditions are imposed at r = Rin by applying the time-dependent analytical solutions. We 
solve the whole region covering 6 £ [0, tt], although our analytical solution are symmetric a,t 6 = Tr/2. So we can check whether 
our code maintains this symmetry. The parameters are the initial time to — 11.0, Rs{t — to) — 10.8, Rc{t = to) — 8.83, 
rg = GM = 0.148, = 1, Ro = 2.58 x 10"^, Oi = 0.1, n ^ I, fi ^ 0, and piRf/A^ = 100. The initial maximum Lorentz 
factor (7s(t — to,r — Rs{t — to))) is 5. Radius and time are normalized by the inner radius Rin and its hght crossing time, 
respectively. 

To obtain the analytical solutions, the ambient plasma pressure po is not specified because we adopted Ranking-Hugoniot 
relations for a strong shock, so that the ambient pressure is negligible. However, in order for the ambient plasma to be in 
hydrostatic equilibrium, the gas pressure gradient force should balance with the gravitational force of the central star. Since it 
is a hard task to reconstruct the self-similar solutions taking into account the ambient plasma pressure, we use the analytical 
solutions obtained by using approximate Rankin-Hugoniot relations given in equations (|48|) - (|50| ) as the initial conditions 
and the ambient pressure is taken initially to be constant, which means that the ambient plasma is not in hydrostatic 
equilibrium. The ambient plasma slowly falls toward the central star due to the gravity of the central star. The parameters 
we used in numerical simulations are taken so that the free fall time t// is much larger than the dynamical time td (typically, 
iff ltd — 10'^). Thus th e free fall motion of the ambient plasma does not affect the dynamics of the expanding magnetic 
loops. Istone &: Norman ( 19921 ) adopted a different method for this problem such that the gravitational force is artificially 
su btracted in t he ambient plasma to numerically recover the solutions of the non-relativistic coronal mass ejection obtained 
bv I Low! (|l984al ). We confirmed that the results obtained by the method proposed by them are consistent qualitatively and 
quantitatively with those including the gravity in the ambient plasma. The ambient plasma pressure, which is not specified 
in the analytical solutions, is taken so small (po = 10~*) that it does not affect the dynamics of the outflows. 

Fig. [6] shows the pressure profile (colour) and the magnetic field lines (curves) a,t t = 48. The left figure shows analytical 
solutions and the right one does numerical results. Fig. [7] shows the density profile (colour). The magnetic field lines in Fig. [6] 
and Fig. [7] are depicted as the isocontours of the fiux function. The levels of the isocontours are identical in both figures. 

As time goes on, the magnetic loops containing the fiux ropes expand in radial direction. The fiux ropes carry the 
toroidal magnetic fields. The ambient matter infiowing through the shock is compressed and accumulated between the contact 
discontinuity and the shock. The numerical results excellently recover the analytical solutions. 
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Figure 7. Colour contour shows the density profile and curves show the magnetic field lines for analytical solutions (left) and numerical 
results (right). 
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Figure 8. The Lorcntz factor (left) and the density and pressure profiles (right) on the equatorial plane at t = 48. Thick curves denote 
the numerical results, while thin curves do analytical solutions. 



Fig. [8] shows the shock structure on the equatorial plane at t = 48. Left panel shows the Lorentz factor and the right 
one does the density and the pressure. Thin solid curves denote the analytical solutions, while thick solid curves do the 
numerical results. At this time, the maximum Lorentz factor is 7 = 7.2 for the analytical solution. The peak Lorentz factor 
in the numerical simulation is, however, 7 = 6.2. The difference comes from the shock flattening in the simulations. Note that 
we adopt the harmonic mean to evaluate primitive variables on the cell surface. This method is more diffusive than other 
interpolation methods, such as the MUSCL type interpolations. Although we utilize this method to avoid numerical oscillations 
at the strong shocks, it decreases the peak Lorentz factor. Also the density and the pressure profiles are diffused (right panel 
of Fig. [8]). Another reason comes from the assumption of strong shocks adopted to derive the approximate Rankin-Hugoniot 
relations (|48p - (|50p . These relations are correct with an accuracy of O^l/j'^). Since we take the initial Lorentz factor at the 
shock as 7 = 5, a few percent error arises from the approximations. 

Such a numerical diffusion produces the sound waves from the shocks. Fig. [9] shows the radial profile of the Lorentz factor 
on the equatorial plane from the initial state at t — to = 11 to t — 25 with the interval St = 2. Thick and thin solid curves 
show numerical results and analytical solutions, respectively, while the dashed line denotes a reference radius r» = 0.995i?s 
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Figure 10. Time evolution of the Lorentz factor at r = r* = 
0.995-Rs. Horizontal axis shows the time and vertical axis shows 
the Lorentz factor on the equatorial plane. Solid curve shows the 
analytical solutions. Grey contours show the Lorentz factors at the 
grid points closest to the reference radius r, = 0.995-Rs. 



Figure 9. Time evolution of the Lorentz factor on the equatorial 
plane. Thick curves show the numerical results, and the thin curves 
do the analytical solutions. Dashed curve shows the reference radius 
r = r. = 0.995/Js. 
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Figure 11. The dependence of the equatorial Lorentz factor at f = 48 on models of the gravity, (i) —GMjry^/r^ (dashed) (ii) 
—GMph'y'^/r^ (dash-dotted) and (iii) without gravity. Solid curve denotes the analytical self-similar solutions. An inset shows solu- 
tions in the different range 1 ^5 r ^ 25. 



(see Fig. I10[). After the simulation goes on, the shock is flattened due to the numerical diffusion, generating the sound waves 
propagating inward and outward from the shock. The compressional waves extract a part of the fluid kinetic energy from the 
shock, resulting in the increase in the numerical velocity behind/ahead of the shock from the analytical solutions. As the wave 
propagates away from the shock front, the radial profiles approach those of analytical solutions. Note that the amplitudes of 
the inward propagating wave decreases with time because the wave conserves the wave energy density = P7^. Since the 
density increases inward, the velocity deviation ^7 decreases as the wave propagates. 

Fig. llOl shows the time evolution of the Lorentz factor. Solid curve shows the analytical solution. We measure the Lorentz 
factor at a reference radius r» = 0.995J?s(t) to avoid the effects of the shock flattening. Since the simulation is carried out at 
discrete grid points, we plot the range of the Lorentz factors of the mesh points closest to the reference radius r, = 0.995i?s 
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Figure 12. The dependence of the equatorial Lorentz factor (left) and the density profile (right) at t = 48 on the density parameter of 
the ambient gas 5. Dashed curves denote the numerical results for initial condition given by the analytical solutions (5 ~ 3.5). The other 
curves depict those for 5 = 5, 4, 3, 2, 1, 0. 



(grey contours). The numerical results deviate from the analytical solution when t < 13 because the flattening of the shock 
front and the emission of the sound waves temporarily increase the velocity behind the shock. The numerical results are, 
however, consistent with the analytical solutions after the sound waves propagate away (t > 18). Then the Lorentz factor 
increases with time as 7 oc t^^*. 

Next we numerically verify the ap proximation of neglecting the thermal enthalp y term /i]v in the gravitational force 



IMext we numerically veruy tne ap proximation or neglectmg tne thermal entnalp y term Aiiv m the gravitational force 
GMph-y'^/r'^ = -GMp{l + hN)y'^/r^ (jlVIobarrv fc LovelaceJllQSd : iMeliani et allbood ). When we derive analytical solutions 



of the self-similar expansion, we neglect in the gravity. Although the specific thermal enthalpy Hm becomes larger than 
unity behind the strong shock (ps/pa — Fa 2> 1, see equations I48I50|) . the gravitational force itself becomes small compared to 
the other forces in the self-similar stage when Rs ^ rg. To evaluate the contribution of the gravity, we carried out numerical 
simulations for three models of the gravity, i.e., (i) —GMp^^/r^, (ii) —GMph'y^/r^ (iii) without gravity. 

Figure [11] shows the radial profiles of the Lorentz factor on the equatorial plane at t = 48 with different models of the 
gravity. Solid curve shows the self-similar solutions, while dashed, dash-dotted, and dotted curves show the numerical results 
for models (i)-(iii), respectively. An inset shows solutions in the different range of r, 1 ^ r ^ 25. Behind the shock, numerical 
results are almost independent of the models of the gravitational force, indicating that the gravitational force is much smaller 
than the other forces. As we mentioned in § 2, the ratio of the gravity for the thermal enthalpy to the plasma inertia decreases 
with radius (see equation I25|). The ratio rg/Rs is 0.01 at the initial state t = to and 3 x 10~^ at the final state t = 48 in our 
simulations. Thus the gravity for the thermal enthalpy is negligible. This explains why the numerical results are independent 
of the gravity models. The specific thermal enthalpy ftjv is larger than unity just behind the shock, but the gravitational force 
is much smaller than the other forces when rg <^ Rs- 

Although the gravity becomes important in the region where Vg < r <ti Rs, the thermal enthalpy is negligible in this 
region (i.e., Hm «C 1). Thus, we can neglect the contribution of /ijv in the gravity (see the inset of Fig. Ill[) . We note that 
the numerical results for model (iii) deviate from the analytical solutions in this region. The plasma is accelerated in radial 
direction by the pressure gradient force, leading to the formation of shocks (r ~ 18.5). Since we use analytical solutions for 
inner boundary conditions at r = Rin — 1, the plasma is supplied from the inner boundary. When the plasma is not confined 
by gravity, the outflowing plasma forms second shocks at r ~ 11. 

When Rs is close to rg, the contribution of /ijv in the gravity is not negligible, but we have to take into account the 
general relativistic effects in such region, so that the characteristic length rg enters into the formulations. In such regime, no 
self-similar solutions can be obtained. It is out of the scope of this paper to obtain solutions in this regime. 

Next we carried out simulations with different density profiles of the ambient plasma to study the generality of the 
analytical solutions and the effects of the ambient density distribution on the loop dynamics. We substitute the density profile 
given in equation (|56l) with the power law profile as. 



-s 

r 



p(r) = po(r = Rs,o) l^j^J , (72) 

where Rs,o — Rs{t ~ to). The analytical solutions correspond to 5 ~ 3.5. We study the several case (5 = 0, 1, 2, 3, 4, 5). The 
initial condition is given by the analytical solutions inside the shock. 
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Figure 13. Time evolution of the Lorcntz factor at the reference radius r*. Circles depict the numerical results for the initial conditions 
given by the analytical solution for 5 ~ 3.5. Other symbols show the numerical results for 5 = 5,4,3,2,1,0. The Lorentz factor is 
evaluated at r» = 0.995r(7,„ax), where r{7niax) is the radius where the Lorentz factor is its maximum. 



The dependence of the equatorial Lorentz factor (left) and the density profile (right) at t = 48 on the density parameter 
of the ambient gas 5 are plotted in Fig. 1121 Dashed curve denotes the numerical results for initial conditions given by the 
analytical solutions (5 ~ 3.5). Other curves depict those for (5 = 5, 4, 3, 2, 1, 0. The peak Lorentz factor decreases as 5 decreases. 
This is because the shell becomes massive for a smaller 5 by sweeping up the larger ambient plasma. The swept up plasma 
is accumulated behind the shock surface. The inertia from the excess plasma accumulated behind the shock decelerates the 
outflows and creates another discontinuity behind the shock. The discontinuity can be considered as the reverse shock. The 
compression ratio of the reverse shock is larger for the denser ambient plasma. 

Time evolution of the equatorial Lorentz factor at the reference radius r« are shown in Fig. 1131 Circles depict the numerical 
results for the initial conditions given by the analytical solutions for 5 ~ 3.5. Other symbols show the numerical results for 
5 = 5, 4, 3, 2, 1, 0. The Lorentz factor is evaluated at r* = 0.995r(7max), where r(7inax) is the radius where the Lorentz factor 
is its maximum. 

The peak Lorentz factor increases with time when t < 13. This increase comes from the emission of the sound waves 
propagating from the shock front. When t > 18, the peak Lorentz factor increases with time when S > 3. The critical value 
of S whether the outflow is accelerated or not can be evaluated from the mass conservation. The plasma density of the 
outflow decreases with radius by according to the mass conservations (see equation I16|l . The rest mass energy of the 
ambient plasma accumulated in the shell thus increases with time when S < 3. On the other hand, when 5 > 3, the outflow 
is accelerated since the plasma inertia of the outflow decreases with time. Naively, we can understand these processes from 
equation (|71[) . According to this equation, the flow is accelerated when the ambient plasma density decreases faster than r~^. 
When (5 = 5, the shock is accelerated and its Lorentz factor is proportion al to the rad i us (~ time). It indicates that the flow 
expands freely. The influence of the ambient plasma is almost neglig ible. IPiran et al.1 (|l993l ) derived the self-similar solution 
of the free expansion. The shock Lorent z factor then increa ses with radius linearly. The shock profile for 5 = 5 or a larger S is 
consistent with the solution obtained by Piran et al. ( 19931 ). Inside the shock surface, our solutions are, however, not identical 
with their analytical solutions since our solutions include the magnetic fields and are intrinsically non-spherical. 



6 SUMMARY & DISCUSSIONS 



We derived axisymmetric relativistic self-similar solutions of the magnetic flux rope expansion by assuming the purely radial 
flow and ignoring the stellar rotation. By taking the self-similar variable as r; = r/Z{t), the arbitrary function Z{t) has a 

unique form given in equation (|22|) . The MHD equations are then solved analytically. 

The solutions obtained in this paper are the extension of our previous work I Takahashi et al. 20091 ) by considering the 
two discontinuities, the contact discontinuity and the shock. The contact discontinuity separates the outflowing plasma and 
the ambient plasma. The outgoing waves propagating in the ambient plasma fo rm sh o cks. 



Such a self-similar solutions including two discontinuities are derived by iLowl (|l984al ) in non-relativistic MHD. Our 
solutions are the extension of their solutions to the relativistic MHD. For the non-relativistic case, the compression ratio 
at the shock is determined by the specific heat ratio. The specific heat ratio is taken as F = 4/3 in non-relativistic MHD 
equations. The system is marginally stable for the inertial flow. The compression ratio is then up to 7 for the strong shocks. In 
relativistic plasma, the sound speed is limited to ~ 0.58 for the ideal gas. The differences between the upstream flow velocity 
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and the downstream wave velocity are larger for the larger flow velocity. This fact results in forming the strong shocks and 
the compression ratio can be larger than 7. The ambient plasma is abruptly heated by the strong shocks. The hot plasma 
is accumulated in the shell between the shocks and the contact discontinuity. Inside the co ntact discont inuity, the magnetic 
loops anchored to the central star are assumed to follow the flux rope solutions obtained bv iLowl (jl984al ). The flux ropes are 
contained inside the global m agnetic loops. S uch magnetic field configuration can be expected for the SGR fiares after the 
magnetic energy is dissipated ( Lvutikov 20061 ) . 

We also carried out numerical simulations of two dimensional relativistic MHD by using the self-similar solutions as the 
initial and inner boundary conditions. Since analytical solutions are obtained in this paper, we can apply them to check the 
accuracy of multi-dimensional relativistic MHD codes. Many previous authors reported that the relativistic MHD code is 
verified by using the one dimensional shock tube problems. There are only a few standard multi-dimensional problems, such 
as the blast wave problem or the rotor problem. However, no exact solutions are known for these problems. We have shown 
that the self-similar solutions can be applied to check the accuracy of the relativistic MHD codes. 

Numerical calculations show that the shock velocity strongly depends on the ambient plasma density. When the density 
profile is steeper than oc r~^, the shock Lorentz factor increases with radius. On the other ha nd, it decreas es for profiles 
shallower than r"''. Such a behavior is expected from the consideration of the mass conservation l|Shapirolll979l ). The density 
inside the shocks approximately decreases with radius as oc r~^. When the decrease in the ambient plasma density is steeper 
than the shocked plasma density, the outflows can be accelerated. We can understand this from the energy conservation given 
in equation (|7ip . The time dependence of the shock Lorentz factor is related to the ambient plasma density. When Rs — t 



and the ambient plasma density is represented by the power law on r (p oc 



the shock Lorentz factor is expressed as 



r, oc r(3-*'/^ 
with radius. 



Numerical results agree with this relation. Especially when 5^5, the shock Lorentz factor linearly increases 



Finally let us apply our results to the magnetar flares ( Woods fc Thompson 20061 : Mereghetti 20081 ). In the discussion 
below, we concentrate on the extraordinarily energetic outbursts (giant flares) observed in SGRs. Although AXPs as well as 
SGRs would be magnetars, the giant fiares have not been detected from the AXPs. The reason would be that the giant fiares 
are very rare events (1 per a few decades in the whole sky). 



The shock Lorentz factor Fs is estimated from the mass ejected by the fiare Mejc 



47rpo^s/3, and equation 



24 



lO-iB erg J \ 1022 g 



(73) 



The expansion speed of the magnetic loop is relativistic when A/cjo 
behind the forward shock is estimated from equations (|49p . (|50p and 



< 10^ 
(1731) Hi 



g (see, also Lvutikov 2006h . The mass density 



Ps — 1.6 X 10 g cm 



f £' 



\ 10« erg 



1022 g 



V109 cm) 



(74) 



The mass density in the ambient plasma po ~ ps/i'^^^'^F s) estimated from equations (|49|). (|50|l . (|73|l and (|74|l . is much 
larger than the Goldreich- Julian density. Such dense coronal pair plasm as would be created by the pair production when the 
magnetically trapped fireball is formed on the surface of the magnetar (jBeloborodov fc Thompsonll2007l ). 
The temperature of the shocked gas is evaluated by using equation (|48l) as 



r, ~ 2.8 MeV 



1046 erg 



M 



1022 



(75) 



Here we assume the pair plasma. Although the radiation spectrum of the initial spike in the giant fiare is not well det ermined 
because of its short duration, the typical temperature indicated by the spectrum is a few 100 keV (jHurlev et al.ll2005l ). which 
is lower than Ts estimated from our analytical solutions. We have to note that the temperature at the photosphere Tps should 
be smaller than Ts because the temperature decreases with decreasing the radius behind the shock (see. Fig. [5] and Fig. [8]). 
The radius of the photosphere _Rps can be determined by the condition that the optical depth of the expanding magnetic loops 



Tps satisfies 



pft7(l 



)dr = 1, 



(76) 



where 9 is the angle between the velocity vector and the direction of the photon propagation (jAbramowicz et al.l Il99ll ) 



and k — ^ ks{kcb + nff) is the effective opacity. Here Kcb and are the opacity for the electron scattering and the free- 
free absorption, respectively. We assume that the ambient plasma (r > Ra) is optically thin. The electron scattering is the 
dominant source for the opacity inside the shocks. This is because the plasma temperature is increased by the shock heating, 
so that the free-free opacity {pT~^'^) is much smaller than that of the electron scattering. 

We numerically integrate equation (|76p assuming that 6 = 0, the expansion energy £' — 10*® erg, and the mass of the 
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central star M = IMq, where Mq is the solar mass. The temperature at the photosphere Tps calculated at R - 
analytical solutions (equations 1 13|[T5] 1161 1531 and I54 p can roughly be fitted by 

when 10 < Fs < 100 and 10^ cm < -Rs < 10^" cm. The result is almost independent of the radius of the contact discontinuity 
(i.e., Rc/Rs)- The temperature in the observer's frame is Tobs = 77ps — FsTps — lOOkeV when Ts — 10. Here we assume 



10 



(77) 



7 : 



Vs since the photosphere is very close to the shock surface. This result is consistent with the observational results when 



10 and Rs 



l(f cm. 



Such hot, relativistically expanding magn etic loops are expected to be formed with the help of the magnetic reconnections. 
According to the scenario bv lLvutikovj ( 2006 ), the ma gnetic reconnections inside the magnetic loops are responsible for the 
initial spike of the flares (see, also Gill fc Hevl 2OI0I ). Subsequently, the plasma is heated up by the shocks produced by 
the magnetic reconnection. By applying our models to the giant flares, the time evolution of the luminosity from expanding 
plasma is expressed as L oc T^\^^R^^ oc V~^'^ R~^'^ from Tobs ~ TsTps and equation (|77|) . Here we assumed R^s — Rs- Since 
Fs oc R^s~^^^'^ from equations (I71|l and (|72|l . we obtain L oc ^-0.8(1+*) {-.gg^ygg Ji^ oi t from equation (|63|l . Terasawa et al 



( 20051 ) reported that the observed photon counts decreased exponentially with time after the initial spike. However, the photon 
counts at the earlier stage, whose decay time is very short (< 100 ms), is not inconsistent with the po wer-law decay. 

In addition to the initial spikes, a hump is observed a few hundred seconds after the initial spike (jTerasawa et al.ll2005l ). 
It is considered that the energy is re-injected from the central star. We now consider another possibility for this hump. Before 
the magnetic energy release, the toroidal magnetic energy can be comparable to that of the poloidal magnetic fields inside the 
magnetic loops. When the magnetic reconnection takes place inside the magnetic loops, the magnetic energy of the poloidal 
magnetic fields is converted to the plasma energies, generating the twisted flux ropes. The toroidal magnetic fleld energy of 
the flux ropes does not dissipate in this process. When the twisted flux ropes cross the Alfven radius (~ light cylinder), the 
rest magnetic energy will be dissipated by interaction with the ambient global magnetic fields. The time scale that a fiux rope 
crosses the light cylinder and releases the magnetic energy is about ~ 1 sec, which is consistent with the observations. The 
energy dissipation at the light cylinder can be responsible for the humps. We need further study to verify these processes. 
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